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ABSTRACT 



C r-v r ». ov LfgnARy 

X'.'O CCMO^L 



The single channel autoregressive lattice has been 
successfully applied to problems including speech analysis 
and recognition, spectral analysis and noise cancelling. 

More recently the two channel autoregressive (AR) lattice 
has been exploited for autoregressive moving average (ARMA) 
analysis of systems for modeling and identification. 

This dissertation considers the multichannel AR lattice 
when applied to ARMA systems analysis. Constraints on 
lattice parameters, based on the input output relations of 
the system under test, are developed. The lattice is 
redefined in terms of the frequency domain representation of 
the input data. This proves to be useful because it allows 
the input to be normalized so that the lattice yields a 
consistant set of parameters independent of the test source 
characteristics. Lastly the lattice is redefined in terms 
of correlations of the input signals . This results in a 
computationally and storage efficient lattice algorithm. 
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I . INTRODUCTION 



In many applications the use of mathematical models 
provides critical insight into dynamic system behavior and 
is essential to many control, simulation and communication 
applications. Models have been applied in the areas of 
economics, neurophysics and astronomy with varying degrees 
of success. Engineers routinely use models when they wish 
to predict or control a system^ performance. 

In this dissertation one form of model, the multichannel 
lattice, is examined in detail in an effort to gain insight 
into its operation and to develop more efficient algorithms 
for its implementation and application to the modeling of 
multichannel linear and nonlinear autoregressive moving 
average (ARMA) systems, based upon a least mean square error 
approach. 

A. THE PROBLEM STATEMENT 

The output y(k) of the discrete, linear, time invariant 
system of Figure 1.1 may be described as the weighted 
summation of the present and past values of the input u(k) 
and past values of the output y(k) . This can be summarized 
by the relation 



q 

y(k) = I a u(k-n) + 
n=0 n 



P 

I 3 y(k-n) 
, n 
n = l 



(1-la) 



3 



defining the transfer function 



Y(z) . q(z) 
U. ( z ) STz) 



where 



q _ n P _ n 

a(z) = E a z and 3 (z) = 1- E 3 z (1-lb) 

n n „ _ n 



u(k) 


Linear Time 
Invariant Digital 


y(k) 




System 





Figure 1.1 General Discrete System 

where a^, CKn_<q and b , lj: n J;P are Fhe parameters of the 
system and k denotes the sample index. 

An estimate of the system parameters can be made if the 
output signal, y(k), and possibly the input signal, u(k) , of 
a system can be observed. The result is a model of the 
system. There are numerous reasons why this might be done, 
among which are: 

1) The parameters of the system are not known and some 
insight into the system function is desired. 
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2) The nominal values of the system parameters are 

known and we wish to detect a change in the system 
performance . 

Equation (1-1) suggests that a system output is predictable 
from a linear combination of past outputs and inputs. 

Because of this, procedures to find these parameters are 
often referred to as linear prediction. 

One way to construct a model is to hypothesize a system 
of the same form, namely 

y(k) = Z a u(k-n) + Z b y(k-n) (1-2) 

_ n n -) n*^ 

n=0 n=l 

A /\ 

and adjust the model parameters a^ and b n so that the 

A 

difference between y(k) and y(k) is minimized according to some 
criteria. Figure 1.2 shows this method. The symbol is 

used to indicate an estimated value. From Figure (1.2) and 
(1-1) and (1-2) 

e(k) = y(k) - y(k) 

q a p k s*. 

- £ (a -a )u(k-n) + Z 3 y(k-n) -Zb y(k-n) 

n-0 n=l n=l 

A 

Since y(k) is itself dependent on model parameters, and 
a minimum mean square output error criterion results in a 
set of nonlinear equations, the solution of this equation 

A 

tor model parameters is generally intractable when b , 

n 

l£n<p is not zero. 
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Figure 1.2 Direct System Modeling 

Because of this, other analysis procedures have been 
developed. A very powerful model is the equation error 
formulation, which entails two models; the analysis model, 
with parameters a^ and which are estimates of the system 
parameters, and a synthesis model using a^ and b n to emulate 
the system. 

The output of the equation error analysis model is 
written as 

q p 

y(k) = Z a u(k-n) + Z b y(k-n) (1-3) 

n=0 n n=l n 



li 




where the estimated output is now a function of the system 
past and present input and system past outputs (instead of 
the model outputs). This form is shown as Figure 1.3. With 
the equation error approach a MMSE solution for the model 
parameters is a set of simultaneous linear equations with a 
unique minimum. The synthesis model for the equation error 
approach is of the form of (1-2), using the estimated system 
parameters found from the analysis model. The equation 
error of Figure 1.3 is given by 

- q p 

e(k) - y(k) - y(k) = E (a -a )u(k-n) + E (3 -b )y(k-n) 

n=0 n n n=l n n 

(1-4) 

which is linear in the model parameters. When a =a and 

n n 

3 n =b n , e(k)=0. In terms of the z-transform, ( 1-4 ) can be 
written as 

E(z) = -A(z)U(z) + B(z)Y(z) 

= C-A(z) + |[|y a( z ) ] U(z) (1-5) 

where 

q p 

A(z) = E a z n and B(z) = 1- E b z. n . 

n n -i n 
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Figure 1.3 Equation Error Analysis Model 

Three cases are examined in this dissertation. 

1) The all pole, or autoregressive (AR) , model 

a = 0 , 1 <n<q 
n — — * 

2) The all zero, or moving average (MA) , model 

b = 0 , 1 <n <p 
n — — 

3) The pole-zero, or autoregressive moving average 
(ARMA) , model 
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a t 0 , for some n, l<n<q 
n — — 

b t 0, for some n, l<n<t 
n — — * 

The objective will be to develop procedures and algo- 
rithms that allow more efficient calculation of these models 
for application to system performance evaluation and 
identification. 

B. OVERVIEW OF THE DISSERTATION 

The following chapters are first a review of some 
contemporary thoughts on linear systems modeling using 
autoregressive, moving average and autoregressive moving 
average modeling, with their lattice formulation, followed 
by some new lattice configurations centered upon the ARMA 
lattice configurations. 

Yule [19] first developed the autoregressive model in a 
paper on sunspot analysis. Because of its wide ranging 
applications, autoregressive modeling has received extensive 
examination since. Chapter II. A begins by reviewing the 
relevant theory of the least mean square solution of the AR 
equations. This solution results in a set of simultaneous 
equations that are the equivalent of the Yule-Walker [9] 
equations. Their resolution requires the solution of the 
matrix equation Rb = r, where R is a symmetric and Toeplitz 
correlation matrix. A number of methods have been proposed 
for the solution of these equations. Matrix inversion 
methods can be used, but their computational complexity is of 
concern for higher order problems. Adaptive schemes [18] 
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have been developed that update an estimate of the solution 
as each data point is taken into account. Levinson [6] 
presented an elegant and efficient recursive method of 
performing this inversion by finding an (n+l)-th order 
solution as an adjustment of the n-th order solution. Durbin 
[2] later recognized that the elements of the vector r were 
also elements of the matrix R. This was used to improve 
the computational efficiency of Levinson's algorithm. This 
recursive solution can be implemented as a lattice structure 
that has been found to offer an accurate and efficient 
solution to the AR modeling problem. More recently Burg 
[1] has shown that this lattice is a maximum entropy 
solution, proving that the lattice solution is optimum. 

In Section II. B the moving average model is similarly 
discussed. The result of a least mean square solution is 
the digital form of the Weiner-Hopf [17] equations and also 
requires the inversion of a symmetric and Toeplitz matrix. 
Perry [16] has shown that this can be efficiently solved 
using Levinson's recursion and can also be restructured to 
form a lattice implementation, part of which is the AR model. 

Sections II. A and II. B deal with single channel (single 
input, single output) models. Section II. C discusses 
multiple input, multiple output generalizations of the AR 
and MA equations and structures that Robinson [16] has 
shown can be solved as a multichannel generalization of the 
single channel structures. 
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Chapter III is a discussion of ARMA modeling. Unlike 
the AR and MA models, the ARMA model represents both poles 
and zeros and so will offer a lower order solution than 
either the AR or MA models for zero-pole systems. In III. A 
the least mean square solution of the ARMA problem is 
developed. A direct application of the ARMA equations 
results in a set of nonlinear, multimodal equations. This 
problem was greatly simplified when the equation error 
formulation [5, 12] was found to result in linear equations. 
Then in Section III.B it is shown how Perry [16] recon- 
figured, with appropriate assumptions, the ARMA problem 
so that it could be solved using a structure identical to 
the two channel AR lattice. This structure offers the 
advantages of a recursive in order solution and results 
in more accurate solutions than batch methods that use 
more direct matrix inversions . 

Lattices have, in the past, used data in the form of a 
discrete time series as a block of data, as discussed in 
Chapter II and III. A and III.B. Griffiths [3] and later 
Morf [13] developed adaptive algorithms that still use 
discrete time series data but update an estimate of the 
lattice parameters as each data point enters the lattice 
structure. In Section III.C a new formulation for deter- 
mining the lattice parameters is presented. This approach, 
although fundamental, has generally not been exploited, 
particularly for the multichannel case . In this structure 
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the data presented to the lattice is the frequency domain 
spectrum of the signal to be analysed. The power spectral 
desnity of the error spectrum is minimized and the resultant 
lattice is equivalent to that used with time domain data. 

This new lattice structure is useful because it allows the 
spectram content of the input signals to be observed and 
presents the possibility of adjusting the input spectra 
to compensate for a nonwhite driving source. This configuration 
also permits the lattice to be used for filter synthesis . 

Section III.D discusses the lattice ARMA model when applied 
to system identification. The ARMA lattice identifier can be 
used for system performance evaluation. In this application 
the model parameters obtained as a result of a test are com- 
pared with model parameters obtained from a system functioning 
normally and functioning with known faults. This comparison 
can reveal if the system under test is operating satisfactorily 
and possible identify a class of faults within the system. 

The lattice parameters fully identify, in a least mean square 
sense, both the input and output of the system to be modeled. 

For modeling it is generally desirable that the input signal 
be white (have a flat spectrum) so that the signal energy is 
equally distributed across the spectrum of interest. Some new 
relations among the lattice parameters are developed, which, 
with some previously known relationships , reduce the number 
of parameters needed by one half when the input is properly 



17 



conditioned. Some new methods are then derived that allow 
preemphasis or normalization of the system input to correct 
for a nonwhite input signal. This is important because 
lattice parameters, while modeling the system, are still a 
function of the input signals second order statistics. This 
is useful even when the test source is purportedly white, 
since the lattice must function with estimates of signal 
parameters obtained from finite length data signals. 

Section III.F introduces another application of the new 
frequency domain lattice. The lattice can be used to experi- 
mentally synthesize a filter to meet a desired specification 
by providing the desired response as the input of the 
lattice. The model generated is a minimum mean square error 
approximation of the desired specification. 

In Chapter IV and Appendix A a new approach to both the 
single channel and two channel lattice is developed. This 
method recognizes that the error signals of the lattice are 
only linear combinations of the input signals and delayed 
versions of the input signals. From this observation the 
equations for the lattice parameters are restructured so that 
only linear combinations of the input correlation functions 
are needed to solve the lattice. This is more computational 
and storage efficient than other lattice methods since it is 
not now necessary to generate updated error signals at each 
stage. It retains the advantages of the lattice of orthogonality 
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and maximum entropy. It offers the new advantage that data 
need not be artifically windowed to force the correlation 
matrix to be Toeplitz, as is required by Levinson's 
algorithm. This may lead, with some types of data, to a 
more accurate solution. This method is extendable to higher 
dimensioned lattices as are used for nonlinear and multiple 
input multiple output system modeling. 

Chapter V presents a summary of the new results in this 
dissertation. Included is a discussion of the application 
and implementation of these new concepts, along with questions 
for future research. 
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II. LINEAR SYSTEM MODELING AND IDENTIFICATION 



The autoregressive, moving average, and autoregressive 
moving average models have all been successfully used for 
system modeling and analysis. The AR model is the model of 
choice when only system output is available. The MA model 
can give lower MSE if an input signal is available, but 
cannot generate an infinite impulse response. These models 
will be examined in greater detail and computationally 
efficient algorithms for their solution will be developed. 
Both single input and single output and multiple channel 
implementations are discussed. 

A. THE AUTOREGRESSIVE MODEL 

For the autoregressive model the assumption is made that 
the present output of a system can be estimated from a 
weighted summation of the past values of the output. This can 
be expressed as 
P 

y(k) = I b y(k-n) (2-1) 

n = l n 

or in matrix form 

y(k) = b T y (k) (2-2) 
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with 



b T = [b n b 0 . . .b ] 

— 12 p 

y T (k) = [y (k-1) y(k-2) ... y(k-p)] 

If a signal y(k) is to be modeled the prediction error 
can be written as the difference in the signal and its 
estimate or 

e(k) = y(k) - y(k) = y(k) - b^y(k) (2-3) 




Figure 2.1 An AR Prediction Error Model 
Figure 2.1 illustrates this type of system. 
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The mean squared prediction error as a function of the 



weights b is found to be 
° n 



E = b x R b - 2 b r + R (0) 

— yy- — yy yy 



(2-4) 



where in general R (n) = e {x ( k) w( k+n ) } , r - e (x(k)w(k+n)} , 

xw — xw — 

T 

R = e{x(k)w (k)} and e{*} designates expectation. 

-xw — — 

In this case, 



R 

-yy 



T 

— yy 



R 


(0) 


. R „(-p+l) 


yy 




yy 


R* 


(p-i) 


. . . ’ R (0) 


Lyy 




yy J 


CR „ 


(i) . . 


. R (d) ] 


yy 


yy * 



The solution of (2-4) for the optimum weights is given by 



b . = R r 
-opt -yy -yy 

which are commonly called the normal equations . 
substituted into (2-4) then the minimum error is 



(2-5) 

If this is 
found to be 



E . = R ( 0 ) - b . r 

mm yy —opt — yy 



( 2 - 6 ) 



(2-3) may be written in the z domain as 



E(z) 

Y(z) 



1 



D 

2 

n = l 



b z 
n 



-n 



3(z) 



(2-7) 



which is the relationship of error signal output to the 
input signal. 
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Consider the all pole filter of Figure 2.2 with the 
input signal u(k) and output y(k) . The relationship 
between input and output can be written in the z domain 



Y(z) 

OTzT 



H( z) 



0 



p 

1- z 

n=l 



= B(z) 



b z 
n 



-n 



( 2 - 8a ) 



or in the time domain 



P 

y(k) = a_u(k) + Z b y(k-n) (2-8b) 

0 , n 




If the output y(k) of the all pole filter is applied to the 

input of the prediction error model of Figure 2.1 the output 

error will be e(k) = a.u(k) . The transfer function of this 

u 

cascade arrangement from the input of the all pole filter to 
the error output of the predictor will be a„ . This shows that 

U 

if u(k) is white noise the effect of the prediction error 
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model is, except for a gain term a , to reverse the action 
of the all pole filter or to whiten it's output. For this 
reason the prediction error model is often referred to as a 
whitening filter or an inverse filter, and makes it useful 
as a spectrum analyzer. It should be noted that a^ establishes 
a lower bound on the MSE of the AR predictor. If u(k) is a 
unit impulse, u(k) = 5 (k), the model may replicate the 
response y(k) perfectly, but can never predict the arrival 
of the input a^5 (k) since the predictor is causal, and there 
is always this minimal error. 

The relationship between an all pole filter of Figure 
2.2 and the analysis model of (2-1) suggests that the all 
pole filter can be used as a synthesis model for generating 
the signal y(k) using a Q u(k) as the input signal with the 
coefficients obtained from the prediction error analysis 
model. It is only necessary to find a suitable gain term 
a^. Equation (2-3) can be rewritten 

y(k)=e(k)+b T y(k) (2-9) 

Comparing this with (2-8b) it is seen that a Q u(k) = e(k) . 

By requiring that the total energy of the output from the 
linear predictor be equal to the total energy of the 
input to the all pole filter, a^ is specified by 



2 _ s {e z (k)} 

0 " R (0) 
uu 



( 2 - 10 ) 
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This leads to a z domain transfer function representation 
for the synthesis model of 

H(z> = (2-11) 

The impulse response of (2-11) can be of infinite dura- 
tion. This means that a relatively low order AR model may 
emulate an infinite impulse response system to an 
acceptable accuracy. 

In the problem as stated so far it has been tacitly 

assumed that the desired order for the model was somehow 

known. In practice this is often not the case. Instead, 

(2-5) is solved for some order and the MSE determined. If 

this MSE is unacceptable, then (2-5) is again solved for a 

higher order. Solution of the simultaneous equations of 

3 2 

(2-5) requires, by Gaussian elimination, p /3 + 0(p ) 
operations, so that a repetitive solution of the normal 
equations can become computationally ponderous. There is, 
however, a high degree a symmetry in these equations and 
the exploitation of this symmetry can give a more efficient 
solution of (2-5). 

1 . A Recursive In Order Solution 

The normal equations (2-5) for the optimum weights 
can be written in matrix form 
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R (0) R (-1) 

yy yy 


. . . R (-p+1) 

yy 




1 

■ — 1 
XI 

1 




R (1) 
yy 


R (1) R (0) 

yy yy 

• • 

• # 






b 2 




V 2) 


R (p-1) . . . 

l_yy 


R (0) 

yy 




1 

Cb 

rO 

1 




1 

*d 

1 



The autocorrelation matrix is Toeplitz (the elements along 
any diagonal are equal) if the signals are stationary which 
is a basic assumption of this work. 

Calculation of the autocorrelation functions now 
deserves some comment. Computation of the exact discrete 
correlation requires a summation of infinite duration 
signals . 

QO 

R (N) = Z y(k)y(k+N) 

yy k =— 



The result is an even function, or 



R 



yy 



(N) 




( -N ) 



so that the autocorrelation matrix is symmetric and 
Toeplitz. In the absence of an infinite duration signal, 
the more common case, an estimate of the correlation func- 
tion must be made. If a window function is used and the 
data is taken to be zero outside the window then the 
estimated correlation matrix will be symmetric when the data 
is averaged over the width of the window. If the data is 
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! 






not taken to be zero outside the window and averaging is 
done over a finite number of points, the estimated 
correlation matrix is in general not symmetric unless a 
large number of data points are taken because of end 
effects. In the derivations which follow symmetric Toeplitz 
matrices are assumed. It is noted that the final algorithms 
can be used without assuming symmetric Toeplitz correlation 
matrices, and apparently work just as well. In fact it is 
noted that Griffith's time adaptive technique [3] is 
nonstationary and yet appears to work very well. This has 
not been explained in the literature. If the observed signal 
is ergodic an estimation of the correlation function can be 
made using a finite interval of data. 

When the data is windowed so that the autocorrelation 
matrix is Toeplitz and symmetric, the method used to solve the 
resultant normal equations is often referred to in the 
literature [9] as an autocorrelation method of solution. 

When the data is not so windowed and the matrix is not 
Toeplitz, the problem solution is said to be a covariance 
method. Either method results in an estimate of the true 
values of the correlation matrix, and it is generally 
unclear that one method will be more accurate than the other. 

The Levinson algorithm relates the (n+l)-th order 
solution of a set of equations of the form of (2-12) to the 
n-th order solution. The assumption is made that b^ n+ '^, 
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where the superscript denotes the order of the solution, can 
be found from b (n ^ by adding a correction to each term and 
one new term, or 



(n+1) _ 



b <n) 


+ 


£. ( n ) 


v , (n) 

where b = 


, (n) 
b i 


0 




k (ntl) 




, ' (n) 

D 

_ n 



(2-13) 



To find b^ n+ ^"^ parameters and must be computed, 

Permuted versions of the vectors b^ n ^ and r^ n ^ may be 

-yy 

defined by reversing the order of their elements . 



(n) 



[>> 




R <n) l 


n 




yy 


• 


(n) 




• 


P = 






-yy 




b (n) 




R (1) 


_1 J 




Lyy J 



= r 



(n) 

-yy 



(2-14) 



The symmetry of the autocorrelation matrix allows (2--12) with 

p=n namely R b = r , to be rewritten using (2-14) as 

- J — yy — — yy ° 



0 (n)~(n) _ (n) 

R b - p 

-yy - -yy 



(2-15) 



The matrix R^ n+ '*'^ is recognized as the matrix R^ n ^ with one 

-yy -yy 

new row and column and may be written as 



R 



(n+1) 

-yy 





R (-n)~ 


' 


r~ 




R (n) 


| : yy 




R (n) | 


(n) 


| R (-1) 




-yy | 


p 


-yy 


yy 

L _ 




T 


— yy 


1 

k < 

( < 

✓ — \ 
2 

k < 


r 

| R (0) 

yy . 




T | 

P (n) 1 

-yy 1 


R (0) 

yy J 
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since the data has been properly windowed to ensure that the 
autocorrelation function estimate is an even function of 
time. The foregoing and (2-13) may be substituted into an 
(n+l)-th order version of (2-12) to give 



R (n) 

-yy 


1 (n) 

-yy 

1 




. (n) , ^(n) 
b +^_ 




r r (n > 

-yy 


T 

p <n) 

-yy 


i 

1 R (0) 

1 yy 




k (n+1) 




R (n+1) 

yy 



Multiplying the matrices and substituting (2-15) yields 



»(n) _ _^(n)] < (n + l) 



(2-16a) 



where 



k (n+l) 



R( n+1) 

yy 

R( 0 ) - 

yy 



- p 

- - y 



( n) 1 , ( n) 



(n) 



T 



Pn) 

-yy 



Substituting (2-16a) into (2-13) yields 



, (n+1) 
b = 


k (nf 


+ k <ntl> 


-b (n) " 




0 




1 



( 2-16b ) 



(2-17) 



m jr • j , (n+1) , (n) , . 

To find b from b it is only necessary to 

determine k^ n '^, which can be found from (2-16b) which uses 

R^(n+1) and other terms from the previous solution. It is 

noted by comparison with (2-6) that the denominator of (2-17) 
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is the MMSE of the (n)-th order solution for the normal 
equations. Equations (2-13), (2-14), (2-16) and (2-17) can 
be applied recursively, starting with order 1, until a model 
results with the desired accuracy. 

2 . An AR Lattice Structure 

The vector b^ n ^ can be written in the transform 

domain 



B (n) (z) = 1 - 



n 

Z 

i = l 



, (n) -l 
b . z 

l 



= 1 - 



h (n) -1 
b i z 



, (n) -n 
-b z 
n 



(2-18) 



and the reversed in order version 




can be written 



3 (n) (z) = 



z -n 3 (n) (z -1 ) 



(2-19) 



It should be recalled that B(z) is the transfer function of 
the prediction error model. The overbar here is used to 
indicate the reversal in order, and will be used throughout 
this dissertation to signify a backward signal, or a function 
or martix associated with a backward signal. 

A recursive solution for the transfer function B(z) 
can be found by substituting (2-17) into an (n+1) order 
version of (2-18) 



B (n+1) (z) = 1 - tz'V 2 .. .z' n - 1 ]b (ntl) 
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(n) 



, r -1 -2 -n-l-i 
1 - L z z . . . z J 



+ k 



( n+1 ) 



-b 

1 



( n) 



B (n) (z) - z- 1 k (n+1) z- n [z n z n - 1 ...l] 






and substituting (2-19) into the second term of the foregoing 
The result is 



B (n+1) (z) = B (n) (z) - z' 1 k <n+1) 3 (n) (z) 



( 2 - 20 ) 



B(n + D(z) can also be found recursively by substituting 
(2-20) into (2-19) rewritten for order (n+1) 



B (n+1) (z) = z 1 B (n) (z) - k (n+1) 3 (n) (z) 



( 2 - 21 ) 



Since B^ n ^(z) is the transfer function of the n-th 
order prediction error model, when it is driven by the 
signal Y(z) the result is the n-th order error signal, or 



E (n) (z) = B (n) (z)Y(z) 



( 2 - 22 ) 



The forward prediction error is the difference in the signal 
and a predicted signal derived from a weighted summation of 
the past n samples of the signal. The backward prediction 
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error is the difference in the signal at time (k-n) and 
its predicted value based on a weighted summation of the n 
samples to follow. Figure 2.3 is an illustration of this 
forward and backward prediction. 




Figure 2.3 Illustration of Forward and Backward 
Prediction Using n Points 



By multiplying (2-20) and (2-21) by Y(z), 
E(z) may also be found recursively. 



E(z) and 



(n+1) 



=r(n+l 



<n) (z> - k ( n+l ) z - 1 E (n) (z) 


(2-23) 


- 1 E (n) (z) - k <n+1) E (n) ( z) 


( 2-23a) 
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In the time domain these two equations become 



e (n+1) (k) = e (n) (k) - k (n+1) e (n) (k-1) (2-24) 

— (n+i) (k) = e (n) (k-l) - k (n+1) e (n) (k) (2-25) 

A 

Since the zero order prediction of a signal is y(k) = 0 (no 
prediction at all) the zero order error is 

e (0) (k) = i (0) (k) = y (k) (2-26) 

and equations (2-24), (2-25) and (2-26) can be implemented as 
a lattice structure as is shown in Figure 2.4. 




Figure 2.4 Second Order AR Lattice Model 
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The lattice parameter k^ ri+ '^ can be found from (2-16b) 
(as was done for the Levinson algorithm) or can be found by 
minimizing the expected value of the squared error of (2-24) 
or (2-25). Minimizing square of the the expected value of 
(2-24) with respect to k yields 



k 



(n+1) 

F 



s (e (n) (k)e (n) (k-l)} 
£{i (n) (k-l)i (n) (k-l)} 



(2-27) 



and similarly for (2-25) 



k 



(n+1) 

B 



e {e (n) Ck>i Cn) <k-l>} 

e (e (n) (k)e (n, (k)} 



(2-28 ) 



These alternative solutions are referred to as the 
forward and backward methods and can be shown [8] to be 
equivalent to (2-16b). From (2-20) and (2-21) it is seen 
that the forward and backward prediction transfer functions 
are equal so (2-27) and (2-28) must be equivalent. However, 
when implemented, the expectations will be estimated from 
time averages and will generally not be equal. This has 
led to two alternative solutions, the first a harmonic 
mean of (2-20) and (2-21) due to Burg [1] 
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(2-29) 




e {e (n) (k)e (n) (k-1)) 



BURG 




and another from Itakaura and Saito [4] which is the geometric 
mean of (2-27) and (2-28) 



Stability of the result is ensured if all the poles of the 
transform domain transfer function lie within the unit 
circle. Markel and Gray [11] have shown that the necessary 
and sufficient condition for this is 



Since (2-30) is of the form of a normalized correlation it 
is of necessity always bounded by unity, and since the 
magnitude of a geometric mean is an upper bound of the 
magnitude of a harmonic mean, (2-30) and (2-29) are 
guaranteed to result in stable solutions . There is no such 
guarantee with the forward (2-27) and backward (2-28) methods. 



(n+1) e (e (n) (k)e (n) (k-1)} 



(2-30) 




|k (n) | < 1 



(2-31) 
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A most important property of the lattice structure 



is the orthogonality of the backward error signals, 
described by 



e (e^^(k)e < '^(k)} 




(2-32) 



0, i*j 



Indeed, the lattice structure can be derived [11] starting 
with this property. This orthogonality decouples successive 
stages of the lattice, ensuring that the addition of stages 
in the recursive solution of the filter will not affect 
preceding stages, and permitting the derivation of (2-27) 
and (2-28). Furthermore, [8] orthogonality allows for a 
recursive calculation of the error energy at successive 
stages by 




p (n+1 ) _ p (n) (1 _, < (n+l) 2 ) 



(2-33) 



Noting that 




(2-34) 



(2-33) and (2-34) may be used for the denominator of either 
(2-16b) , (2-27) or (2-28) . 
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Once (2-17) is solved, the vector b^ n ^ may be 
implemented as the all pole filter of Figure 2.2, or other 
equivalent arrangements, to produce a synthesis model. 

The lattice synthesis structure may be generated by (2-25) and 
rearranging (2-24). 

e (n) (k) = e (n+1) (k) + k (n+1) e (n) (k-l) (2-35) 

This equation can be realized by the structure of Figure 2.5 
and implements the transfer function 1/B^ n ^(z). 




Figure 2.5 Second Order AR Lattice Synthesis Structure 



B. THE MOVING AVERAGE MODEL 

The moving average model is based on the assumption that 
the present value of a system output can be estimated by a 
weighted summation of the present and past q values of the 
system input. In equation form this is 
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(2-36) 



y(k) 



q 

Z a(n) u(k-n) = a 
n=0 



T 



u(k) 



where 



a T = [a(0) a(l) ... a(q)] 



u~(k) = Cu(k) u(k-l) ... u(k-c)] 



(2-37) 



It should be noted that the vectors of (2-37) are (q+1) 
elements long. 

Defining the error of the model as the difference in the 
estimated value and the actual value of the system output, 
or 



where R is (q+l)x(q+l) and r is a (q+1) element vector. 

— uu n n — uy 

When this is minimized with respect to the coefficients of a 
the result is the normal equations of the moving average 
problem 



e(k) = y(k) - y(k) 



(2-38) 



the mean squared error E is found to be 



E = a R a 
uu — 



T 




(2-39) 




(2-40) 
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From this the optimum value for a is 



a 4 . 
—opt 




(2-41) 



while the MMSE is found by substituting (2-41) into (2-39) 



E . 
mm 



R (0) 



yy 



T 

3 . 

—opt — uy 



(2-42) 




Figure 2.6 Moving Average Modeling 

Figure 2.6 shows the block diagram of the implementation 
of the MA model. The transfer function of the model is, in 
the z domain 

Q -1 

A(z)=Za(n)z (2-43) 

n=0 
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and is often referred to as an all zero model. As such its 
impulse response is finite and in all cases the model is 
stable . 



1 . A Recursive In Order Solution 

Equation (2-40) is similar in structure to (2-12). 

R is Toeplitz as is R of (2-12), but the elements of r 
— uu r — yy — uy 

of (2-40) are not also elements of the matrix R . Still 

— uu 



Levinson's algorithm may be applied. 

Equation (2-40) is written for order n as 



D (n) (n) (n) 

R a = r 
— uu — — uy 



( 2-44) 



and the assumption is made that a^ n+ '^ is a function of a^ n ^ 
and a correction vector. 



(n+1) 



£ <nf 


+ 


Y ( n+ D 


0 




s <n + l> 

o 



(2-45) 



Define vectors 



r (n+1) = 

— UU 


hr (i) 

uu 


(n+1) 

, P = 

— uu 


Hr ( n+ 1 r 

uu 




R (n+1) 
uu 




R (1) 
_ uu 



~ (n+1) 

— UU 



(2-46) 
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and expand an (n+1) order version of (2-44) with (2-45) 
substituted . 



-< n ) 
— uu 



(n+1) 

p 

— uu 



(n + 1 ) 

p 

— uu 



R (0) 
uu 



(n) (n+1) 

a + y 



g 



(n+1) 





(n) 

-VO 




-uy 




R( n+1 ) 




L u y A 



When this is multiplied two equations result 



(n) (n) + (n) (n+1) + (n+1) (n+1) = ^(n) 

-uu - -uu - -uu 6 -uy 



(2-47) 



(n+1) (n) , (n+1) (n+1) , p rn s (n+1) _ . 

p a + p y + R (0) g - R (n + 1) 

— uu — — uu — uu 6 uy 



(2-48) 



The first and last terms of (2-47) are equal so both may be 
cancelled leaving, when solved for Y^ n+ ^ 



Y 



(n+1) 



-R (n) 
— uu 



-1 



(n+1) (n+1) 

P 2! 

-uu 



(2-49) 



Finding R^ n ^ Y^^ is comparable to (2-5), the (n + l)-st order 
autoregression of the input signal u(k) . So to derive the 
MA model, the AR model of the input must first be solved. 

When this is done a new vector f^ n+ D can ^ e defined where 
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f ( n+1 ) 



(n) ~ 
— uu 



(n+1) 




(n+1 
n+ 



i 



b 



(n+1 ) 
1 



(2-50) 



is the result of the autocorrelation analysis, 
defined as 



Now 



Y (n+1) 



is 



(n+1) p (n+l) (n+1) 

1 = - £ g 



(2-51) 



and this can be substituted into (2-48) and solved for 



(n+1) 

g 



g (n+1) 



^(n+l) (n+1) (n) 

R - p a 

uy — uu ^ 

R (0) - p (n+1)T f (n + 1) - 
uu — uu — 



(2-52) 



2 . A MA Lattice Structure 

The Levinson solution of (2-49) can be restructured 
to give a lattice configuration. The transfer function of 
the MA model (2-43) can be written recursively in the z 
domain using (2-45) and (2-51) as 

A (n+1) (z) = A (n) (z) + g (n+1) B (n+1) (z) (2-53) 



4 2 



where B^ n+ ‘*'\z) is the result of the autoregressive analysis 
of the input u(k) . If both sides of (2-53) are multiplied 
by the input U(z) and then transformed to the time domain the 
result is the estimate of the output signal y(k) for 



~ (n+1) 

y 



(k) 



V n) (k) + g (n+D— (n+l) 



(k) 



(2-54) 



where e^ n+ ^(k) is the backward error signal of the autore- 
gressive analysis of u(k). The desired lattice structure 
is shown in Figure 2.7. 



e 



( 0 ) 



u(k) 




Figure 2.7 A Second Order MA Lattice Structure 
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Recalling that the backward error signals of different 
stages of an AR lattice are orthogonal, it is seen that the 
MA estimate of y^ n ^(k) is a linear weighted sum of these 
orthogonal signals. 

C. MULTICHANNEL MODELING 

The models previously considered have been of single 
input, single output systems. Multiple input, multiple 
output systems can be modeled as vector generalizations of 
the models heretofore developed. 

1 . A Multichannel Generalization 

A single channel output of a Q-channel autoregressive 
model is the weighted summation of the past outputs of the Q 
channels of the system to be modeled. Similarly, a single 
channel output of a Q-channel moving average model of a 
system is the weighted summation of the past N inputs . of -.the 
Q channels to the system to be modeled. Either of these 
statements may be written in equation form 

Q N 

x.(k) = E E d..(n) x.(k-n) (2-55) 

3 i=l n=l 13 1 



Q = number of channels 
n = time delay index 
i = input channel 



k = time index 

j = output channel being 
estimated 
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A 

Xj(k) = estimate of output channel j at time k 
x^(k-n) = input to channel i at time (k-n) 

d.-(n) = weighting factor for estimating the output signal 
1 -' at channel j at time k due to input x (k-n) 

As an example, for a model with 3 inputs (Q = 3) , 

(2-55) would be expanded as 



N 



N 



N 



x^(k) = Z d-^(n) x ]_^ c-r 0 + Z d 2 ^(n) x^(k-n)+ Z d 3 ^(n) x^(k-n) 
n=l n=l n=l 



N 



N 



N 



x 2 (k) = Z d^ 2 (n) x^(k-n)+ Z d 22 (n) x 2 (k-n)+ Z d 32 (n) x 3 (k-n) 
n=l n=l ^ n=l 



N N N 

x^(k) = Z d^ 3 ^ n ) x^(k-n)+ Z d 2 ^(n) x 2 (k-n)+ Z d^ Q (n) x 0 (k-n) 
n=l n=l n=l 



Equation (2-55) can also be written in matrix form. 

x'(k;) T = d~ X (k) = [X (k) T D] T (2-56) 

where 

x(k) T = [x^(k) x 2 (k) ... x^(k)] (Qxl) 

X^k) = [x^(k) x‘(k) ... x£(k)] (NQxl ) 

xT(k) = [x i (k-l) x ± (k-2) ... x ± (k-N)] (Nxl) 
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ill 


-12 


• • • q in 


—21 

f 


i 22 


* * • d 2N 


1 — 1 

. -JT 


^Q2 


t 

‘ • ' d QN 



(NxNQ) 



di . = [d . . ( 1) d . . ( 2 ) ... d . . (N) ] (Nxl) 

-13 13 13 13 



The error of such an estimate is the vector difference 
in the actual output vector and the predicted output vector. 



e (k ) = x ( k ) - x(k) = [I - D 1 ] X(k) 



The prediction error covariance matrix is a (QxQ) matrix 



P = e[e(k) e (k) ] 



( 2 - 57 ) 



If the trace of the prediction error covariance matrix is 
minimized the result is 

RD = r ( 2 - 58 ) 
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where 



R = £ (X (k) X (k)} = £ 



x x ( k )~ 



x^(k) 



[ x ^ ( k ) 



• • Xq(k) ] 



R 



“- 1-1 



R 

^Q-l 



• 5. T 

— 1— Q 



. R 



( 2 - 59 ) 



is a matrix of autocorrelation and cross correlation matrices 
and 



_r = £ 


x^(k) 


[x^(k) ... Xp(k) ] = 


r 

-1 X 1 


... V 

' — 1 X Q 




x r ,(k) 




^Q X 1 


... r 

X Q X Q 



( 2 - 60 ) 

" is a cross correlation matrix of size (NQxN) and elements 



r = e 

x . x . 

-i ] 


X. (k-1 ) 
i 


x . (k) = 
] 


R / x 

x . x . ( 1 ) 
i ] 




• 




x i (k-N) 




R (N) 

X .x . 

i : 
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As imposing as this may seem, it is no more than a 
vector generalization of (2-12) from the AR problem and 
(2-40) from the MA problem. Each correlation coefficient is 
replaced with a correlation matrix in the multichannel 
generalization . 

It will be useful in the sequel if the two channel 
lattice is examined additionally. Equation (2-56) expanded 
for two channels is 



x 1 (k) 


T 


i — 1 
i — 1 

| Ttl 


' H 

CM 

T3l 




x^k) 


A 


= D X (k) = 










x 2 (k) 




_~12 


CM * 
CM 




x 2 (k) 



N 

= £ 


d,, (n) 


d 21 (n) 




x (k-n) 


2 

II 


c 

v-/ 

CM 

i — 1 


d 22 (n) 




x 2 (k-n) 



The two channel error is the vector difference in the signal 
to be estimated and the estimation. 



e(k) = 



1 — 1 

1 * 




1 ^ 

M 

1 — 1 

|< X 




A 

(x^(k)-x (k) ) 


x 2 (k) 




.x 2 (k) 




(x 2 (k)-x 2 (k) ) 



( 2-63) 
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Equation (2-58), the result of a minimization of the error 
equations, is in two channels. 



R 

X 1 X 1 


* OJ 

X 
1 — 1 
X 
& 1 




-ii 


in' 




r 

- x i x i 


-X X X 2 P 


Csl 

X 

i — 1 

X 

p^l 1 


R 

X 2 X 2_ 




i 2 i 


-2 2 




r 

-x 2 x 1 n 


r 

- x 2 x 2 J 



( 2-64) 



where 



R 

— x . x . 
i ] 



R (0) . . . R (1-N) 
x i X j X i X j 



R (1-N) . . . R (0) 

X • X . X . X . 

1 3 1 ] 





R (-1) 

x i x i 




R (1) 

X 1 X 2 


r = 

x l x 2 


R ( -N ) 

X -x . 
l l 


E Vip = 


R (N) 

X 1 2 



X 1 X 2N 



R 



X 1 X 2 



R ( -N ) 
X 1 X 2 
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Important to the solution of (2-58) is the form of 

the (NQxNQ) correlation matrix R. The matrices on the main 

T 

diagonal are all correlation matrices of the form R and 

— 1—1 

so are symmetric and Toeplitz. The off diagonal 
matrices are cross correlation matrices that are Toeplitz 
but not symmetric. Symmetrically opposed matrices are, 
however, transposes of one another, with the result that the 
matrix R is also symmetric. This block Toeplitz form is 
amenable to a block form of Levinson algorithm solution. 

The vector generalization of the AR lattice is shown 
in Figure 2.8 and is described by the equations 

e (n+ ! ) (k) = e (n) (k) - K (n+1) e (n) (k-l) (2-65) 



i (n+1) (k) = e (n) (k-l) - K (n+1) 



e (n) (k) 



( 2 - 66 ) 



and the initial conditions 

e (0) (k) = e (0) (k) = x(k) (2-67) 

The signal paths of the single channel AR model have been 
replaced with parallel vector signal paths, adders with 
vector adders and multipliers with matrix multipliers. 
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Figure 2.8 A Two Stage Multichannel AR Lattice 

Lattice predictor coefficients and may be 

also found from vector generalizations of the single channel 
case using 



K 



(n+1) _ p(n) a/ 1 "*) 



( 2 - 68 ) 



^•(n+1) _ p(n) 1 ^ (n) 



(2-69) 



A (n) = e (e (n) (k) e (n) (k-l)} 



(2-70) 



,(n) _ _r (n),. . ( n) 



P ww = e (e ww (k) e Vi1 ' (k)} = P* 11 ” 1 ^! - K (n) K (n) ] 



(2-71) 



P (n) = e{e (n) (k) e (n) (k)} = F (n “ 1) [i - K (n) K (n) ] 



(2-72) 
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As in the single channel case, the optimum predictor 
for all orders less than N are found with the N-th order 
solution. Also, orthogonality applies to the backward 
channel error signals, or 



_.T 

e {e 1 ( k ) e • - 1 ( k ) } 







i = j 



and successive stages are thereby decoupled. 

For a multichannel MA model equations (2-65) to (2-72) 
are used with a vector generalization of (2-54). 






t G (n+i)i e (n+1) (k) 



(2-73) 



The error of the MA model is expressed by 

e^ n \k) = y(k) - y (k) 

— o — — 

and can be minimized with rhe result 

Q (n) _ p(n) -1 e{e (n) (k) £ T (k)} (2-74) 

Figure 2.9 illustrates the resultant structure. 

In the next chapter the multichannel generalization 
is applied to ARMA modeling. 
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III. THE AUTOREGRESSIVE MOVING AVERAGE MODEL 



IN TIME AND FREQUENCY DOMAIN 

While the MA and AR models are suitable for many appli- 
cations, modeling some systems can require a very high 
order model to obtain suitable results. Consider the 
geometric series expansion. 

00 

1 - n -n | — 1 1 

— = 1 — = E a z |az | < 1 

(1 - az ; n=0 

This shows that a single pole is equal to an infinite number 
of zeros, and conversly a single zero is equal to an 
infinite number of poles . Thus a system with a single pole 
will require an AR model of infinite order, as would a 
system with a single zero require an infinite order MA 
model, to be precisely represented. The ARMA model emulates 
both poles and zeros , and so can model pole-zero systems 
without resorting to infinite order solutions. 

This chapter reviews contemporary ARMA lattice struc- 
tures, then develops some new formulations using frequency 
domain data and considers their implementation as a 
performance evaluation tool. 
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A. THE ARMA EQUATIONS 



To model a system with both poles and zeros it is 
intuitively appealing to construct a model that also has 
both poles and zeros. As was earlier discussed, the equation 
error model will result in linear equations as the MMSE 
solution. If the condition that the order of poles equal 
zeros is imposed the equation error model output estimate is 
described by 



y(k) 



N N 

Z a u(k-n) + Z b y(k-n) 
_ n n n' 



( 3 - 1 ) 



[y X (k) | u^(k) ] 



where vectors y(k) , u(k) , b and a are defined and dimensioned 
as in the AR and MA discussions of Sections II. A and II. B 

A 

and y(k) is the ARMA output estimate. This constraint is 
not limiting since some of the coefficients of the numerator 
or demoninator can be zero. Minimizing the mean squared 
error with respect to the parameter vectors a and b yields 
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(3-2) 



with a MMSE of 



T 

'"min 



R (0) 

yy 



T I T 
opt | — opt' 



R (if 

yy 



R (N) 

-Tl— 

R (0) 
uy 



R (N) 

_ u y . 



( 3-3) 



Equation (3-2) can be readily solved for and a^^ by 

matrix manipulation but again a recursive in order solution 
■■ can be derived . 



B. A RECURSIVE IN ORDER ARMA SOLUTION 

The ( ( 2N+l)x( 2N+1) ) correlation matrix of (3-2) is 
neither Toeplitz nor symmetric, bur is nearly so. The 
offending elements are those associated with the a^ term. 
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If it is assumed that a Q can be found by other means, (3-2) 
can be restructured so that the correlation matrix is block 
Toeplitz by subtracting the a Q terms from both sides of the 
matrix equations (3-2). 
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(3-4) 



Since a^ is just the DC gain, and can be easily calculated, 
this is a reasonable assumption. Equation (3-4) can be 
rewritten using the notation of (2-64), as was used in the 
discussion of the two channel AR lattice in Section II. C 
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— yu — uu 

— J n - 1 




Q» ; 



(3-5) 



where the vector a contains the N elements a.,... a,. of the 

— IN 

T 

vector a. Here use has been made of the identity R =R~ . 

— J — uv — vu 

This is compared with (2-64) with x^ = y and x 2 = u 
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To solve (3-5) using the algorithms for the two channel 
lattice it is only necessary to make the substitution 
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1 OvJ 

1 1 

'Ol 
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a 




f 21 


—2 2_ 




-a 

o 



(3-6) 



after solving for the vectors cL ^ in order to find b and a. 
Since a Q is the DC gain of the system to be modeled it can be 
found using 



_ s {y ( k ) u ( k ) } , q 7 

a 0 ~ e (u(k)u(k)} 

1 . An ARMA Lattice Structure 

Since the equation error AR'HA model can be solved 
in terms of a two channel AR recursive formulation, it is 
not difficult to modify the AR two channel lattice to 
represent the ARMA system. 

The ARMA error can be written as the difference 

A 

between the output y(k) and the estimated output y(k) or 
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(3-8) 



e Q (k) 



y ( k ) - [ y * ( k ) i u T ( k ) ] 




from the error equation model of Figure 1.3 and (3-1). The 
AR error can be written by substituting (2-62) into (2-63), 
with y(k) for x^(k) and u(k) for x 2 (k), and transposing, 
so that 



e(k) 



T 




[y(k) T u' (k) T ] 



-11 ^12] 

-21 —2 2 j 



= [e (k) e (k) ] " 

y u 



(3-9) 



where u f (k) is an N element vector containing the elements 
u(k-l) . . .u(k-n) of the vector u(k) . 

If (3-9) is post multiplied by a vector ip_ where 

4 =[\J (3-10 



the result is 



e^(k)ii; = [e (k) e (k)] 
— — y u 




= y(k)-a u(k) - [y(k) u ,T (k)] 




e 0 (k) 



(3-11) 
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which is exactly the ARMA error of (3-8), Furthermore, the 
mean squared ARMA error is 

e{eQ(k)} = s{e(k) e(k) T }^ 
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L 0 J 



= P t (3-12) 

where P is the forward prediction error covariance matrix of 
the two channel AR model. If (3-12) is minimized with 
respect to a^ the result is 

e {e (k)e (k)} 

a n = ^ i (3-13) 

u e (e^(k)J- 

which is an alternative solution for a n . 

To construct an ARMA lattice it is only necessary to 
derive the ARMA error using (3-11) from the AR lattice. 

A second order lattice model is shown as Figure 3.1. 

From (2-65), forward prediction in the two channel AR 
lattice can be written 



e (k) 
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k 21 
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(3-14) 
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From this equation e^ n ^(k) and e^ n ^(k) can be obtained 



e (n) (k) 

y 



e (n+1) Ck) + k n ( ^ +1) e (n) Ck-l)+k^ +1) e (n) (k-1) 
y 11 y 21 u 



(n) ,, s 
e (k) 
u 



(3-15) 

(n+l),,N , , (n+1)— (n) n , >. , v (n+l)— (n) ,, , 

= e (k) + k. 0 e (k-l)+k 00 e (k-1) 

u 12 y 22 u 



( 3-15 ) 



Equations (3-15) and (3-16) describe the synthesis 
structure that can be used to implement the ARMA model as 
shown by Figure 3.2. For an N-th order realization an input 
is needed for e^‘^(k). From (3-11) 

e (N) = e (k) + a e (N) (k) (3-17) 

y 0 0 u 



and if the ARMA model is an accurate representation e^(k) 
will be small so that (3-17) can be rewritten 



e (N) (k) - a„e (N) (k) (3-18) 

y 0 u 

This provides the input to e^(k) shown in Figure 3.2. 

C. LATTICE ALGORITHMS IN THE FREQUENCY DOMAIN 

In the foregoing it has been tacitly assumed that the 
data available for analysis is in the form of a finite 
length time domain signal. This may not be the case. The 
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data could be in the form of the Fourier transform of the 
time domain signal, or it may prove advantageous to 
process the data in that form. Makhoul [9] has shown that 
linear prediction can be done equivalently in either time 
and frequency domain, but there has been in the literature 
no exploitation of this for the lattice algorithms. Thus 
it seems useful to examine the lattice structures when the 
data is represented as a frequency domain spectrum. 

1 . The Single Channel AR Lattice 

The time domain lattice was formulated using the 
time domain equations (2-24) and (2-25). With frequency 
domain data available the single channel AR lattice can be 
formulated using (2-23) and (2-23a). Its structure is shown 
as Figure 3.3. A unit time delay becomes a multiplication 
by with the input signal represented by the continuous 

spectrum YU). Note that in the backward channel the signal 
E U) is a delayed version of E' vn 'U). The frequency 
domain error of a zero order prediction is 



E (0) U) = E ,(a) U) = YU) (3-19) 

From the forward error channel recursion (2-23) the magnitude 
squared error of order (n+1) can be written as 



_(n+l), . r (n+l), \_ r — i(n) . . . (n+l)=r(n), .Trr.ln))' .. 
r. U) E U) = U U)-k E U)JLE U) 



, (n+l)=( n) , x 
-k E U)J 



(3-20) 
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Figure 3.3 A Second Order AR Lattice in the Frequency 

where E* indicates complex conjugate. The total magnitude 
squared error can be found as the integral of (3—20) 

/ |E (n+1) (to) | 2 do) = / | E (n) (u>) | 2 d 

-00 -CO 



. (n+l)r 

-k L 



[E (n) (a))E (n) 



(w)]du+/ [E 



(n) (oj)E (n) (w)]dtoJ 



+ k 



(n+D f °° j-(n) (w) |2 daj 



(3-21) 



Since XZ* = [X*Z]* the bracketed term of the right side of 
(3-21) can be rewritten as 



65 



00 



0Q 



[ / [E (n) (o))E (n) (a))]da) + f [E (n) (co)E (n) (a3)]da)] 



. rr , ( n ) f N„(n) 

- f IE ( a) )E 



(0))]d03 + 



CO 5* <* 

f [E (n) ( 03 )E (n) (a))] d(D 

OO 



(3-22) 



Since Y(o)> is the spectrum of a stable and causal sequence 
then Re{Y(o))} is an even function and Im{Y(a)» is an odd 
function. This is also true for the error signals E^ n \a)) 
and E^ '(a)). The result of the integrations of (3-22), must 
therefore be real and the two terms of (3-22) are equal, or 
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(3-23) 



Substituting (3-23) for the middle term of (3-21) 
and minimizing the result with respect to k^ n+1 ^ yields 
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or equivalently 
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periodic m-point 
spectrum vr ]<_> (3-26) 

1 mi 

with period mT 



If the discrete correlation theorem and Parseval's 
theorem are applied to (3-26) the result is 



(n+1) _ e{e (n) (k)e (n) (k-l)} _ ,, (n+1) 

- j " k f 

e (e (n) (k)} 



(3-27) 



which is equivalent to the forward solution (2-27). 

The same procedure can be applied to the backward 
error recursion (2-23a) with the result 
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(3-28) 
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(3-30) 



which can be shown to be equivalent to (2-28), the backward 
equation . 
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( 3-31) 



2 . The Two Channel ARMA Lattice 

The multichannel lattice can be similarly developed 
using frequency domain transformations of the time domain 
equations. Doing so, the two channel lattice equations 
become 




i 
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( 3 - 33 ) 
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( 3 - 36 ) 
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K (n+1) = pCn) -1 A (n) T 



( 3 - 37 ) 



^(n+1) _ p (n) _1 A (n) 



( 3 - 38 ) 
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( 3 - 40 ) 



where the limits of integration depend upon the periodicity 
of the available data. If the data is in a discrete form 
the integrations are replaced with summations. 



discrete form the result of a correlation obtained for a 
time shift other than zero will differ when calculated in 
the frequency domain from the result normally obtained from 
time domain data. Consider first two data signals y(k) and 
x(k) , both defined for l_<k_<N , and k an integer. The 
correlation of delay t is obtained by the summation 



It is important to note that when the data is in a 
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and has the effect of windowing the data so that data points 
outside the window are ignored. If these data signals are 
transformed to the discrete frequency domain with a digital 
Fourier transform the resultant spectrum is periodic, which 
implies that in the time domain the data signal has also 
been made to be periodic. A correlation obtained in the 
frequency domain is now 

R (&) = i— / Y(o)) X(-co) e^-j^ dw 

yx 2 it N 

J -IT 

which when transformed back to the time domain is equivalent 
to 




(£). 




y (k)x(k+£ ) 



N 

+ Z y(k)x(k+£-N) 
k=N-£+l 



1 

N-Z 



This is because the effect in the time domain of the fre- 
quency domain time delay e-i^- is a circular shift of the 
data in the time domain. As was discussed in the pre- 
ceding chapter time domain data would normally use a linear 
shift and then be either windowed or would be considered of 
infinite length. Either method introduces it's own dis- 
tortion into the estimation of the correlation. The circular 
shift of the frequency domain correlation is a third and 
different kind of distortion. This is important when the 
number of data points is small, because a bias will be 
introduced into the estimate of the correlation. 
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D. SYSTEM IDENTIFICATION WITH THE ARMA LATTICE 

The problem of system identification concerns the 
determination of a mathematical model of a system from its 
input output relationship. An impetus for this research has 
been the issue of fault detection within a system. If a 
model can be constructed of an operating system, either on 
line or under test, the model parameters might reveal if the 
system is functioning normally and if not, what class of 
component is faulty. While the ARMA lattice has charac- 
teristics which make it particularly attractive for this 
application, particularly its robustness, there are some 
characteristics that still deserve attention. 

In system identification the lattice model would be 
applied as shown in Figure 3.4. The lattice parameters , 

K^ n ' ) and a^ that result from a test can be compared with a 
dictionary of parameters obtained under normal and faulty 
conditions. Using lattice parameters directly, though, 
results in more model parameters than if the vectors a 
and b were used for identification. Identification using a 
and b uses (2n+l) parameters for an n-th order system while 
lattice parameters result in (8n+l) parameters. Therefore 
it is beneficial to reduce the number of lattice parameters 
that must be examined in the identification problem. 
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Figure 3.4 ARMA Lattice Used for System Identification 

1 . Restraints on ARMA Lattice Parameters 

ARMA lattice parameters are a function of input 
signal second order statistics. The input signals u(k) and 
y(k) to the ARMA lattice are related by the transfer function 
of the system under test. This suggests that the ARMA 
lattice parameters are interrelated. To examine the pro- 
blem the equations for the lattice parameters will be 
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expanded in terms of the system input u(k) and output y(k). 
The vector defining the initial conditions for the ARMA 
lattice first stage is 

= e (0) (k) = e (0) (k) (3-41) 



X(k) = 



y (k) 
u(k) 



The forward and backward prediction error covariance 
matrix become 
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(3-42) 
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and the matrix A becomes 

A (0) = e[e (0) (k)i (0)i (k-l)] 
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(3-43) 
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Now consider the result when the signals u(k) and 



y(k) are from a system as shown in Figure 3.5 
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Figure 3.5 System Excited with White Noise 



From linear system analysis theory several statements 
can be made about the second order statistics of y(k) and 
u(k) , assuming that u(k) is white noise. 
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These identities can be substituted into (3-42) and 
(3-43). Thus 
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(3-50) 
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With these, the lattice parameters can be found 
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From (3-51) and (3-52) some restraints on the 
lattice parameters are evident. 
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(3-54) 





(3-55) 




(3-56) 



It can thus be seen that it is not necessary to 



calculate all eight lattice parameters. Only three, plus 
the DC gain a , are needed. Alternatively, the parameter 
matrices can be calculated and either (3-54) or (3-55) used 
for an alternative equation for calculation of a^ . 

Equation (3-53) can be interpreted in light of the 
lattice structure. kj^ and k^^ are the coefficients that 
are the predictors of the future value of u(k) . Since u(k) 
is uncorrelated from sample to sample no prediction can be 
made and these must be zero. 



stages. The forward error vector of the first stage is 



Some of these relations can be extended to later 
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because of (3*-5 3) while the backward error vector is 



e (1) (k-1) = 



F 2 {y (k-1 ) , y (k-2 ) , u(k-2)} 
F 3 {y(k-2), u(k-l) , y(k-l)} 



(3-58) 



where F^{-} indicates a linear function of the arguments, 
is now defined as 



( 1 ) 



.( 1 ) 

A 11 

e {u(k) • F^ {y (k-1) ,y(k-2) ,u(k-2)» 
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a 12 
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x2 



(3-59) 



because of the independence of samples of u(k) . When this 
is applied to (2-65) for k “ the result is 
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(3-60) 
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3y induction this argument can be extended to 



arbitrary order with the result 



• ( n ) . ( n) _ — a i 

k ^2 - ^2 2 — n - 0, 1 5 



(3-61) 



when 



R ( Jl) = a 5(3,) 
uu u 



Using (3-61) the forward error vector can be written 
for order n 



( n) 
e 



(k) 



agU(k)+F u (y(k-l) , . . . ,y(k-n) ,u(k-l) , . . . ,u(k-n)} 



u(k) 



(3-62) 



and this can be used to find the forward error covariance 
matrix 



p (n) 
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If this is now used with (3-59) of order n to find 



the backward Dararaeter matrix the result is 



K (n+1) 



= P 
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F 2 ( n f| 
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( 3-64) 



and from this it is seen that 



k<n) 

K 21 



t-( n ) 

- a 0 k ll 



(3-65) 



u( n) - 
K 22 



t-( n) 
~ a o k 12 



(3-66) 



For the single channel AR lattice the forward and 
backward coefficients were equal, but for the multichannel 
lattice it is known that the forward and backward parameter 
matrices are not equal. Nuttal [7] has shown that the 
determinants of the forward and backward error covariance 
matrices are equal , and from this it can be shown that the 
determinants of forward and backward lattice parameter 
matrices are equal. 
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so it follows that 



It is known that P ^ ^ = P^ 

Det(P (0) ) = Det(P C0) ) 

Also 

p(n+l) _ p(n) j-p_T 7 (n+l) ^(n+D-j 

p(n) ^(n+1) |-^(n+l)~‘*' _ ^(n+D-j 

and 

p(n+l) _ p(n) rj K (n+1) ^(n+l)-j 

p(n) |-^(n+l) -1 _ K (n+l) j ^(n+l) 

It is known that 

Det(AB) - Det(A)Det(3) - Det(B)Det(A) 

•' so that the determinants of (3-68) and (3-69) can 

Det(P (n+1) ) = Det(P (n) )Det(K (n+1) )Det(K (n+1) - 

Det(P (n+1) ) = Det(? (n) )Det(K (n+1) )Det (K (n+1) - 



(3-67) 

( 3-68 ) 

( 3-69) 

be written 
( n+1) ^ 

(n+1) ^ 
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which shows that Det(P (n+1) ) = Det(P (n+1) ) if Det ( P ( n) ) =Det ( P ( n) ) . 
In view of (3-67) it can be inductively reasoned that 

Det(P (n) ) = Det(P (n) ) (3-70) 

Now the determinats of the lattice parameters defined by 
(2-68) and (2-69) can be written 

Det(K (n+1) ) = Det(P (n) )Det(A (n) ) 

Det(K (n+1) ) = Det(P (n) )Det(A (n) ) 

T 

and since Det(A) = Det(A ) 

Det(K (n+1) ) = Det(K (n+1) ) (3-71) 



A summary of these restraints is presented as Table I. 
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SUMMARY OF RESTRAINTS ON ARMA 
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2 . Uniqueness of the Lattice Parameters 

A particularly troublesome consideration in using 
ARMA lattice parameters for system identification is the 
lack of uniqueness of the lattice parameters for the system 
under test. This is evident because the lattice parameters 
are not just functions of the system under examination, but 
also of the test signal used to drive the system. Thus for 
different test signals such as an ensemble of finite white 
noise input signals, an ensemble of lattice parameters are 
obtained. Interestingly, experiments have shown that the 
transfer function coefficients calculated using these lattice 
parameters will generally have less dispersion than the 
lattice parameters from which they are calculated. 

Thus, the selection of a test signal influences both the 
accuracy of the model and the expectation of the lattice 
parameters that result. The mean square value of equation 
error of the ARMA model can be expressed in integral form as 




( 3 -7 2a ) 



-t: 



The equation error can be written in the z domain as 



ECz) = [3(z)H(z) - A(z)]U(z) 



( 3 -72b) 
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Substituting (3-72b) into (3-72a) results in 



_ 1 _ 
2 7T 



TT 



/ 

-7T 



[BCe 



3“ T )Hte juT )-A( e j “ T )]| 2 p(e> T 



) dwT 



(3,-7 3) 



Equation (3-73) shows that the power spectrum of the driving 
source is a weighting function for the resultant modeling 
error. To obtain a model that is equally accurate throughout 
the spectrum, a test source must be used which has equal 
energy throughout the spectrum. An obvious source with a 
uniform power spectrum is a white noise source. Besides 
giving equal weighting to the error spectrum, white noise 
involves, as has been shown in the preceeding section, 
simplified calculations and a simplified lattice structure. 
This, however, is not the total answer. When applied to 
system identification, the correlations and cross correlations 
of y(k) and u(k) generally will not be available. Instead 
these quantities will be estimated from finite sample 
functions of y(k) and u(k) . Even though these signals may 
be ergodic, the values obtained are only estimates of the 
corr el at ion s . 

Experiments with a white noise source have shown 
that different ensembles of input signals will result in 
different lattice parameters for each member of the ensemble. 
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In order to improve the ensemble averages obtained for the 
lattice parameters, some form of input preprocessing can be 
performed to whiten the input or equivalently flatten it's 
spectrum. 

Consider first Figure 3.6a, which is a modification 
of the system of Figure 3.4. An adjustable filter with 
transfer function H(w) has been added to the input of the 
system under test. The transfer function H(u>) can be adjusted 
to made U'(u>) white if the spectrum of the noise source is not 
white and has not zeros on the u> axis. If the system under 
test is linear, then the lattice filter derived from Figure 
3.6b is equivalent to the one derived from Figure 3.6c. 

Figure 3.7 is a block diagram for calculating the frequency 
domain ARMA lattice model with four preprocessing techniques. 

a) Preprocessing in the time domain on both the input 
and output signals of the system under test. 

b) Preprocessing in the frequency domain on both the 
input and output signals of the system under test. 
Since time frequency domain operations are equivalent 

they can be interchanged. For comparison they are considered 
separately to determine if there are any computational 
advantages of one over the other. It should be noted that 
as long as equivalent operations are performed on both input 
and output, as discussed earlier, the resultant lattice model 
will still model the system under test, though the specific 
lattice coefficients may differ. 
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Figure 3.6 Lattice Input Preprocessing 
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Lattice Parameters 

Figure 3.7 Frequency Domain ARMA Lattice with Preprocessing 
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A random noise generator can be represented as shown 
in Figure 3.6c. If this generator provides the source signal 
of Figure 3.7, then preprocessing to make U'(w) white is 
equivalent to the multiplication in the frequency domain of 
U(w) by G~^(w) , assuming that the inverse exists and is stable. 

In the discrete n--point frequency domain U(w) can be 
written as a vector 

U T = [U n Ik , . . , U , ] ( 3-74a) 

— 0 1 n- 1 

G(co) in its discrete form can also be written as a vector 

G T = [G n G, . . . G . ] ( 3-74b) 

— 01 n-1 

Multiplying by G can be expressed in the discrete 

frequency domain as a multiplication of the vector U by the 

T T 

square matrix H, namely U’ =U H, where 

H = [h..] and h..=0 for i i j (3-75a) 

— 13 i] 

and 

h ii = ^i-1 for i = 2 ; •••. n (3-75b) 

where U n . is the expected value of the ith harmonic of the 
input signal U(w) taken over an ensemble of the input signals. 
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For H to be useful G(w) should have no poles or zeros 
on the unit circle. G(w) can be known apriori (or estimated) 
from an autoregressive analysis of the signal generator 
output, or by averaging the spectra of several (.possible 
overlapping) segments of the generator output. Equations 
(3-75a) and (3-75b) are useful when dealing with an ensemble 
of input signals . Alternately for a single member of the 
ensemble 

h^ = Ik ^ for i = 1, 2, ..., n (3-75c) 

can be used. This is equivalent in the time domain to 
exciting the system with a unit impulse . 

From Figure 3.6c U(io) = N(oj)G(w). From Figure 3.6a 
the inputs to the lattice filter are given by 

U’(w) = U(w)H(u>) = N(o>)G(u)H(a>) = N(w) 

when 

H(co) = G^^kw) . Also, 

Y'Cco) = U ' (co)T(to) = NCu)T(u) 



and the transfer function derived by the lattice model is 



Y’ (o>) 
TTcoT 



T(u) 
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From Figure 3.6b U ' (<u) = U(u>)H(u>) and Y* (u)=Y(u)H(w)=U(u)T(w)H(u) . 
The transfer function for the lattice is 



Y' (w) 
U' (u) 



T(w) 



which is equivalent to that of Figure 3.6a. In Figure 3.6b 
the spectra presented to the lattice filter are not functions 
of G(w) , and the coefficients obtained can be used directly 
to identify the system under test . 

E. FREQUENCY DOMAIN ARMA LATTICE SIMULATIONS 

To verify some of these concepts a FORTRAN program was 
written that implements the ARMA frequency domain lattice 
with preprocessing introduced as shown in Figure 3.3. A 
Hamming window can be selected for' preprocessing in the time 
domain. The selectable frequency domain processing imple- 
mented is that of (3-75c). The spectrum Y(w) is divided by 
the spectrum U(u>) to give Y’(u)), the new system output 
spectrum. The input spectrum U(m) is replaced with a unit 
magnitude real spectrum. 

When preprocessing is not desired then the inputs to 
the ARMA frequency domain lattice are changed to U*(n)= 

U(n) and Y*(n)=Y(n). 

For comparison a program was written that solved the 
ARMA normal equations by Gaussian elimination. Listings 
of these programs are available from the author. 
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igure 3.8 ARMA Lattice Simulation Program Block 
Diagram with Input Preprocessing 
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More than twenty different systems , of up to fourth 
order, were used as the system under test. The number of 
data points used for a single analysis varied from 32 to 
4096. The number of points was always an integral power of 
two because of an FFT routine used in the frequency domain 
lattice . 

To illustrate the ability of the frequency domain 
lattice to accurately model a system two examples, typical 
of all the problems exercised are shown here. System 1, 



TABLE II 
SYSTEM 1 



Transfer Function Coefficients 



Numerator Denominator 



A(0) = 


0 , 
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A( 1) = 
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characterized in Table II, is a fourth order low pass system. 
Figure 3.9a and 3.9b show the magnitude of the transfer 
function of the model obtained by exciting the system with a 
single (256 point) sample of white noise from the Gaussian 
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Figure 3.9c Gaussian Elimination 256 Points 
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elements and the frequency domain lattice using no pre- 
processing. As has been observed in the past [15] the 
lattice results in a closer model of the system. Figure 3.9c 
and 3.9d show the pole zero plots of the resultant transfer 
functions . 

Figure 3.10a and 3.10b are the transfer function magnitude 
found with the same methods using 4096 data points. Here the 
lattice model is still a more accurate model than the batch 
method . 

While in most cases the lattice offered superior perfor- 
mance to the Gaussian elimination method. System 2, Table III, 
is a counter example that shows that this is not always the 
case. In this example the system is a second order band pass 
filter. Figure 3.11a and 3.11b are the transfer function 
magnitude plots found using 128 points, and Figure 3.11c and 
3. lid show the pole zero plots. Here it is seen that the 
batch method produces a superior model. Figure 3.12 shows 
the same system analysis using 1024 points, and both methods 
accurately model the system. 

To determine the ability of the input preprocessing to 
compensate for nonwhite input to the system the arrangement 
of Figure 3.13 was used to generate a colored input signal. 

The plant used is that of System 1. The test input is a 
white Gaussian noise source followed by a single pole filter, 
with the pole position varied from -0.5 to 0.5 for different 
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Figure 3.10c Gaussian Elimination 4096 Points 
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Figure 3.10d Lattice Model 4096 Points 
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Figure 3.13 Input Preprocessing Test 
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test runs. The filter was operated under three conditions 
of preprocessing; 

1) No preprocessing. 

2) Inputs U(oj) and Y(w) normalized as in (3?-75c) 

3) Multiplication of the input data by a Hamming window 
and then normalizing in the frequency domain as in 

C 3-75c) . 

Figure 3.14 shows two typical plots of the first stage 
lattice parameters. The horizontal axis is the pole position, 
b, of the input filter of Figure 3.13. The vertical axis is 
the value of the lattice parameter, in this case k^^ and 

for cases 1) and 2). For all first stage parameters, 
case 3) results in essentially identical results to case 2). 

It is clear that the preprocessing has made the lattice 
parameters independent of the input filter characteristics. 

Figure 3.15 shows two similar plots for second stage 
parameters k^“ and k^ • Again the preprocessing is 

effective, and case 3) is essentially identical to case 2). 

( 3 ) 

Figure 3.16 shows parameter k^ . The first plot is 
that of case 1) and 2). Here the parameter shows some 
sensitivity to the pole position, though much less than the 
unnormalized input signal. The second plot is for cases 
1) and 3). The windowed and compensated lattice still shows 
no sensitivity to the input filter pole position. Of the 
eight parameters, this one showed the most variation with the 
pole position. 
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( 4 ) 

Figure 3.17 is the fourth stage parameter • Here the 

parameter calculated without windowing varies over an even 
wider range than the parameters calculated without any 
preprocessing. Several other of the parameters of the fourth 
stage do the same . The parameters of the windowed and com- 
pensated lattice, case 3), are all still insensitive to the 
input pole position. 

It is apparent from this experiment that input signal 
preprocessing is effective in compensating for the effects 
of a nonwhite input signal, and results in a consistent 
set of lattice parameters for identification. The pre- 
processing should consist of windowing of the input data 
(to eliminate the Gibb's phenomenon) and normalisation to 
whiten the signal. 

F. FILTER SYNTHESIS WITH THE ARMA FREQUENCY DOMAIN LATTICE 

So far the lattice filters have been developed as models 
for systems with measurable input and output. It is possible 
to define the input signal to the filter as a unit impulse 
so that the lattice models a desired response in the frequency 
or time domain. In this application the recursive in order 
nature of the lattice solution becomes of interest. The 
lattice can be expanded in order until a suitable fit of the 
desired specification is obtained. 
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Consider what happens when the frequency domain ARMA 
lattice is used to model a system when the excitation of that 
system is a unit impulse. The configuration is shown in 
Figure 3.18. U(n), the spectrum of the input to the system 

under test, is a flat spectrum of unit magnitude, exhibiting 
the previously discussed desirable property of equal energy 
across the spectrum. The output spectrum of the system Y(n) 
is exactly H(n), the discrete Fourier transform of the 
impulse response of the system under test. 




Figure 3.18 Frequency Domain Modeling 



ill 



The lattice parameters obtained as a result of this 
analysis could be used to construct a synthesis model. This 
synthesis model would exhibit the best fit, -in a least mean 
square sense, that a filter of its order could obtain. The 
order of the lattice can be increased to obtain an arbitrarily 
small MSE . When designing to a specification it is often 
desired that some portions of the frequency characteristics 
of the resultant filter be weighted more heavily than other. 
This too can be accomplished with the ARMA lattice. U(n) and 
Y(n) can both be multiplied by a weighting function which 
provides preemphasis to those portions of the spectrum for 
which a more accurate fit is desired. 

1 . Filter Synthesis Examples 

To illustrate this capability it was postulated that 
a digital filter was desired with a magnitude 

| H( n) | = 0.495 cos • n) +.505 , 0 <n<1023 

This desired magnitude was transformed to a minimum phase 
spectrum using a digital approximation of the Hilbert trans- 
form [14]. This spectrum was used as the input Y(n) to the 
frequency domain lattice , with a unit real spectrum applied 
as the U(n) input. The magnitude of the first and second 
order models of this desired spectrum are shown in Figure 
3.19. The second order model magnitude appears to match 
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114 



the desired magnitude, and as shown by the pole zero plot, 
Figure 3.19c, is stable and minimum phase. The digital 
transfer function of this second order filter is 



H(z) 



0.302504 + 0. 495032 z _1 + 0.202512z 
1 + 0 . 00004z -1 + 0 . 000005z~ 2 



As another example the frequency domain lattice is 
used to find a digital equivalent of the analog low pass 
circuit of Figure 3.20 which exhibits 40 dB attenuation at 
5 kHz. This circuit was analysed using the SPICE program, 
a circuit analysis program. This program produced an AC 
analysis of the circuit, showing the complex output when the 
circuit is driven by a unit magnitude signal for frequencies 
from DC to 10 kHz, the equivalent of a Bode plot of the 
circuit transfer function. 

This complex spectrum was applied to the frequency 
domain lattice as in the previous example. The third, 
fourth and fifth order solutions are shown in Figure 3.21. 

The transfer function of the fifth order solution is 

-. 00469-0. 0269 z _J --0. 310 z~ 2 -0.0 3164 z“ 3 -0. 316 z -4 +.00096z 
H(z) = 

1-2 .159z _1 +2 .534z' 2 -1.675z" 3 + .6178z' 4 -0 .0979z -5 
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Figure 3.20 Analog Low Pass Filter 
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The pole zero and phase plots of the fifth order 
solution are shown in Figure 3.21. There is an additional 
pole not shown on the plot that is located at 3 = — 4 - 5 . _ The 
phase plot, Figure 3.21e, shows that the analog circuit has 
a phase at the Nyquist frequency of approximately tt/ 4. Of 
course the digital filter cannot match this phase, it must 
be either zero or tt at the folding frequency. It is 
believed that this discontinuity in the phase is the cause 
of the digital filters non minimum phase characteristic. 

The filter is stable, however, and does meet the desired 
magnitude specification. 
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Figure 3.21b Fourth Order Model 
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IV, AUTOREGRESSIVE LATTICE PARAMETERS CALCULATED 



USING THE CORRELATION MATRIX 

In the foregoing chapters several methods of solution of 
the AR and ARMA normal equations have been discussed. These 
equations can be solved by matrix methods but for large 
order problems this is a computationally expensive method. 

If the available data is windowed to ensure a Toeplitz and 
symmetric correlation matrix the equations can be more 
efficiently solved using Levinson's algorithm. From 
Levinson's recursion the lattice structure can be imple- 
mented but requires more calculations because it is necessary 
to update prediction error signals at each stage. Makhoul 
[10] developed an algorithm that implements the single 
channel AR lattice but does not require the generation of 
these error signals, and thus is more efficient than other 
methods . 

In this chapter a new implementation is developed for 
both the single channel and the two channel AR lattice that 
implements the lattice structure without the expense of 
calculating error signals. In the final formulation there 
is no requirement that the correlation matrix be Toeplitz 
or symmetric, removing constraints on how data is windowed 
and allowing solution of the problem when the system under 
analysis is not stationary. However, since the formal proof 
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depends upon the lattice formulation which assumes symmetric 
Toeplitz matrices, two solutions are presented, one 
assuming Toeplit symmetric correlation matrices , the other 
omitting this restriction. This algorithm maintains the 
lattice advantages of maximum entropy and robustness that 
have made the lattice an attractive analysis structure. It 
requires less storage than conventional lattice implementations 
and can be extended efficiently- to multichannel structures as 
are used for nonlinear systems analysis. 

A. THE SINGLE CHANNEL AUTOREGRESSIVE LATTICE 

In the single channel autoregressive lattice structure 
of Figure 4.1 the forward and backward error signals for 
each stage are updated using the equations from Section II. B. 



( n+1 ) 



(k) = e 




(4-1) 



e 




(4-2) 



/ n ) 

The reflection factors k can be found from the backward 



equation (2-28) 
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Figure 4.1 Single Channel Lattice Structure 




present and past values of the input y(k) , the numerator and 
denominator of (4-3) can be expressed as the weighted sum of 
correlations of the input vector y. An algorithm based upon 
this observation can be expected to be computationally more 
efficient than previous algorithms if the correlations of the 
input signal are known since it is not necessary to update 
the error sequences at each lattice stage. 

From Figure 4.1 the error signal e^ n ^(k) can be written 



as 



e^ n \k) = [y (k)y(k-l) . , .y(k-n) ] w^ n ^ 




(4-4) 



w 

n 



where 




(n) _ £ . (i), (i-1) 

= L k k 




i=l 



i = 0 



w 




(4-5) 
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w[ n) is the summation of the gain factors of all 
paths with delay i from the input to e K J (k) . A 
formulation for w^ n ^ is now developed. Equation 
be written 

e <n) <k) = y (n)T (k) w (n) 
where 

T 

(k) = [y(k)y (k-1) . . . y(k-n) ] 



and 



(n) r (n) (n) (n) 

w = Lw w, . . . w J 
o 1 n 



It should be noted that the y^ n ^(k) and w^ n ^ are 



long. 



Similarly e^ n ^(k) can be written 



e (n) (k) 



= [y (k)y(k-l) . . .y (k-n) ] 
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recursive 
(4-4) can 

(4-6) 



(4-7) 



(4-8) 
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or 



— (n) ,, v _ (n) ,, v " (n) 

e (k) = y (k)w 



( n) . 



(n) 



where w ww is the inverted version of w . With 
notation, a delay of the backward error sequence 
sented by a shift downward of the elements of the 
vector w k , with a zero shifted into the top, or 



e (n) (k-l) = y (n+1) (k) 



_ 0 _ 

’ ( n) 



w 



(n) 



(k) or 



The weight vector associated with e 
be found recursively. Substituting (4-6) and (4- 
equations for the signals e v ; (k) and e^ ; (k-l), 
forward error recursion (4-1), yields 
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Comparing (4-12) with (4-6) rewritten for order (n+1), yields 



(n+1) 

w 




k (n +1 ) 




Similarly, substituting (4-6) and (4-11) into (4- 
backward error recursion, results in 
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Comparing this with 



(4-11) gives the recursion 
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It should be noted that the intital conditions for the 
weight vectors and are 

w (0) = [1] w (0) = Cl] (4-16) 



As an illustration, from Figure 4,1 it can be seen that 
the second stage weight vectors can be written as 
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Applying these to (4-13) yields 
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and to (3-15) 
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(4-20) 



yielding the third stage weight vectors. 

The numerator of equation (4-3) can now be rewritten 
using the definition of the error signals in terms of the 
weight vectors 
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(4-20) 
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It is seen that the numerator of (4-3), required for calcu- 
lation of k , can be calculated with no presumption of 
symmetry or requirement that the correlation matrix be 
Toeplitz. The lattice method may still be used if the 
signals under analysis are not stationary or estimates of 
the correlation matrix are available that do not result in a 
Toeplitz and symmetric matrix. However, if the correlation 
matrix of (4-20) is Toeplitz and symmetric the matrix 
multiplication may be written as summations. 
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These summations are recognized as convolutions of the 
weight vector , or the correlation of the forward and 

backward weight vectors, 

Defining correlation functions of w^ n ^ and as 
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$ (n) C-n-1) = 0 



(4-23) 



Equation (4-.20) can now be rewritten as 
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This expression can be written more compactly as the product 
of two vectors, one with elements of R , the other of 
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The denominator of equation (4-3) could be found using 
similar arguments, but generally offers no advantage to 
previously described methods. The recursion 

£ {e (0) (k)e (0) (k)} 

e (e (n) (k)e (n) (k)} 

should be used. Equation (4-3) may now be rewritten using 
(4-20) or (4-26) and (4-29). 
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1 . Implementation of the Single Channel Algorithm 
The lattice can be implemented using correlation 
weights by applying (4-16), (4-13), (4-27), (4-20) or 
(4-28), and (4-29) and (4-30). A flowchart of the procedure 
is presented as Figure 4.2. 



initial conditions of (4-16) is used. For 
later stages (4-13) gives the recursive solution 
for w (n) . 

(3) The numerator of (4-30) is found using the weight 



For an mth order analysis the procedure will be: 
(1) For stationary or windowed signals find the 



correlation vector r^ n \ For nonstationary 
systems or infinite unwindowed data find the 
correlation matrix R^ n \ 



(2) Find w^ n ^. For the first stage of analysis the 



vector and either a correlation vector r^ n ^ or 
the correlation matrix R^ n ^. 
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Figure 4.2 Flow Chart for mth Order Analysis 
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(4) Find ? (rU from (4-28). For n=0, P Co) = R (0); 

yy 

otherwise use the recursion of (4-28). 

C5 ) ^tn+l) is now found from (4-30). 

(6) If higher order analysis is desired, increment 
n and return to ( 2 ) . 

Figure 4.3 gives a comparison of the number of 
multiplications required to calculate reflection factors for 
a seven stage lattice that uses 256 data points, and assumes 
a stationary system with the data windowed to ensure a 
Toeplitz structure. 



CONDITIONS: 7 stages 

256 data points 

Data windowed for Toeplitz and symmetric 
correlation matrix 



TIME DOMAIN LATTICE: 
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B. THE TWO CHANNEL LATTICE 

The two channel lattice developed in II, C can similarly 
be redefined in terms of correlation weights , the details of 
which are presented in Appendix A. The following is a 
summary of the steps and equations necessary to implement 
the algorithm. 

To implement the two channel algorithm it is necessary 
to calculate eight weighting vectors, two for each of the 
four forward and backward error signals. These can be found 
using the initial conditions 
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The vectors of (4-33) and (4-34) each have two 
which are also vectors, and the matrix multiplicat 
executed as a (2x2) matrix times a (2x1) vector. 

The (2x2) matrix a' 
vectors and a correlation matrix R' 
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The four elements of the matrix are defined by the 

relations 
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matrices are found using the recursions 



The P (n) and P (n) 
previously cited, namely 
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The lattice parameters are found using the A v , P v 
~*( n ) 

and P^ J matrices from 
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A flow chart of this implementation is presented as 
Figure 4.4. 
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Figure 4.4 



The m-'order Two Channel Implementation 
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1. Correlation Lattice Simulations 



To verify the algorithms presented in the forgoing a 
FORTRAN program was written. This program implements the 
flowchart of Figure 4.4. The correlation matrices generated 
by this program are both Toeplitz and symmetric. The 
examples here are the same two systems used in Chapter III.E. 
and represent a run which is typical of the ensemble of tests 
that were made. The distribution over the ensemble was minor 
since measurement noise was not considered. The fourth order 
model of System 1 is shown with Gaussian elimination using 
256 points, in Figure 4.5. Figure 4.5 shows the results 
using 4096 points. Figure 4.5 shows the same comparison of 
System 2 using 128 data points, and Figure 4.8 is System 2 
modeled with 1024 points. In almost all cases the results of 
the correlation method are essentially identical to the 
results of the Gaussian elimination method. This should not 
be unexpected. Both matrix manipulation and correlation 
lattice methods use identical data, the auto and cross 
correlations of the input data signals. Both methods window 
the data identically. Both methods are solving the same 
normal equations, though through differing algorithms. 

The number of multiplications saved by the correla- 
tion lattice, when compared to the conventional lattice 
implementation, is approximately 8p, where p is the number 
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of data points. This is because the error signals are not 
multiplied by the lattice parameter matrices and K^ 1 "^ 
in the correlation implementation. 
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Figure 4.5d System 1, 256 Points, Correlation Lattice 
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Figure 4.6d System 1, 4096 Points Correlation Lattice 
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System 2, 128 Points, Gaussian Elimination 
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Figure 4.7b System 2, 128 Points, Correlation Method 
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Figure 4.7d System 2, 128 Points, Batch Method 
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Figure 4.3a System 2, 1024 Points, Gaussian Elimination 
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Figure 4.8b System 2, 1024 Points, Correlation Lattice 
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Figure 4.8c System 2, 1024 Points, Gaussian Elimination 
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Figure 4.8d System 2, 1024 Points, Correlation Lattice 
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V. SUMMARY, CONCLUSIONS AND OPEN QUESTIONS 



The objective of this dissertation has been to develop 
efficient methods of implementing the AR lattice. Two 
particular applications for system identification and 
modeling have provided the motivation for this research. 

The first is system fault identification. Computationally 
efficient algorithms must be implemented if the fault 
analysis is expected to be done on a real time system. 
Adaptive lattice structures require only moderate computation 
for each data point and have often been implemented in real 
time. The block structures, such as the conventional time 
domain or frequency domain lattices, require greater 
numbers of calculations in bursts for each block of data. 
These structures could be implemented using both parallel 
and pipeline procedures: parallel calculation of the 

correlations of each channel and pipelined calculation of 
the stages of the model . 

The correlation lattice structure also encourages 
parallel processing implementations. In this structure 
calculation of the correlation matrix requires the greatest 
computational effort. These correlation estimates can be 
performed in parallel as the data becomes available. 

The second application considered is the modeling of 
nonlinear systems. For this application the foregoing 
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issues are of even greater importance. Hybrid signals, 
products and cross products of the nonlinear system input 
and output, are generated and the resultant signals applied 
as the input of a multichannel lattice. In this case the 
computational and storage savings of the correlation lattice 
become of major significance and allow the implementation of 
more complex models. It has been found that the two channel 
lattice when applied to ARMA modeling can be simplified 
because of the relationships that exist between the input 
and the output of the system under test. The lattice can 
be simplified further when the system driving signal is white. 

Some of the constraints found are given physical inter- 
pretations that relate to the system under analysis. There 
has been no such interpretation for the relationship of 
the lattice parameters by a factor of the DC gain of the 
system being modeled. It is felt that this simple relationship 
should have some physical interpretation related to the 
system under test. 

The lattice was redefined using the frequency domain 
representation of the input data. This was found to be 
useful because operations in the frequency domain allowed 
for normalization of the input signals to whiten their 
spectra. This whitened spectra provides a standard input 
that allows lattice parameters to be used for identification 



153 



and ensures the previously mentioned relationships of the 
lattice parameters . 

There are two normalization procedures proposed. The 
first requires the averaging of ensembles of the inputs to 
obtain an expectation of the energy spectrum of the input. 
The second method uses the spectrum of a single sample 
spectrum as a normalization function. Either of these 
methods will fail if the energy at any single spectral 
harmonic is zero. 

It is not understood how either of these methods affects 
the accuracy of the resultant model. Experiments comparing 
normalization methods, with and without spectral smoothing 
introduced by windowing, have been inconclusive. While some 
examples have shown normalization can improve model accuracy 
over an unnormalized analysis , other examples show the 
opposite. It is postulated that the dynamic range of the 
input spectrum and normalization spectrum are of major 
significance, but there has been no analysis of the 
procedures involved. 

It is also found that this frequency domain lattice can 
be used to synthesize a digital filter. The desired 
frequency domain spectrum is provided as the input of the 
lattice, and recursive in order lattice algorithm provides 
an LMS lattice model which meets the desired specification. 
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Lastly it is found that the lattice can be reformulated 
in terms of correlation weights of the input signals. This 
is efficient because it becomes unnecessary to generate 
error signals within the lattice. This results in both 
computational and storage efficiencies when compared to 
conventional lattice implementations. Though developed for 
the two channel lattice , the procedure can be extended to 
higher dimensioned lattices where the efficiencies will be 
of even greater significance. 

Inherent in the derivation of Levinson's algorithm is 
the requirement that the input signals correlation matrix be 
symmetric Toeplitz, requiring that the signals be stationary 
and correlations estimates are from windowed data. When the 
lattice structure is derived from Levinson's algorithm the 
requirement of windowing the data is displaced by the 
lattice structure, which forces upon the data signals a 
different but nevertheless present window function. The 
symmetric Toeplitz nature of the correlation matrix is not 
an issue with the lattice structure. 

When the lattice is reformulated in terms of correlation 
weights, it is again necessary to generate a correlation 
matrix. It does not appear in the derivation that there is 
a requirement for a symmetric Toeplitz matrix but it is also 
not understood how a nonsymmetric Toeplitz correlation 
matrix will affect the accuracy of the solution. This 
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formulation can be expected to generate errors in the 



solution, and it is not clear how they will propagate 
through the correlation lattice algorithm. This is an 
that deserve future attention. 
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APPENDIX A 



The two channel lattice structure of Figure 2.8 can be 
reconfigured for a more efficient implementation using 
correlation weights rather than error signals. The lattice 
equations earlier defined in II. C are repeated here for 
convenience . 
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y(k) 
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e (0) (k) 




y (k) 


e (0) (k) = 
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_e (0) (k)_ 




_u(k)_ 



( A- 2 ) 



e (n+1) (k) = s <n, 0O - K (n+1) e (n) (k-1) 



CA-3) 



e (n+1) (k) = e (n) Ck-1) - K (n+1) e (n) (k) 



(A-4 ) 



and the reflection factor solutions 



P (3) = p (0) = £{e ( 0 ) (k)e ( 0 )X (k)} 



R (0) R (0) 
yy yu 



R (0) R (0) 
_ uy uu 



( A-5 ) 
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p ( n+1 ) _ p ( n ) ^ j ^ ^ ^ ^ ^ K ^ ^ ^ J 

p(n+l) _ p (n ^[T,K (n) K (n) J 

. (n) , (n) ,, s— (n) -, s ■, 

A_ = e ie (k)e (k-1)} 

K (n+1) _ p ( n ) 1 A (n) T 
^(n+1) _ p(n) -1 A (n) 



( A- 6 ) 



(A-7 ) 



(A-3) 



( A- 9 ) 



( A-10 ) 



A channel error signal e^ n \k) or e^ n \k) is a linear 
combination of the present value and the past n values of y 
and u. Using this, ep ^(k) may be written 
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or 



e <n) (k) = [y (n)T (k)!u (n)T (k)] 



w<"> 

-yy 



(n) 

w 

-yu 



■ (n n 

-uy 



(n) 

w 

— uu 



159 
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where 



T 

y^ n ^ (k) = [y (k) . . . y (k-n) ] (A-13) 



U (n) (k) = [u(k) . . . u(k-n) ] (A-14) 



and 
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— xz 



= [w 
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. w j 

xz 

n 



(A-15) 



Vectors y^ n ^(k) u^ n ^(k) and are each (n+1) terms long. 
The result of the multiplication of (A-12) is a 2 element 
column vector. 
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Using this same notation the backward delayed error 



vector e (k-1) may be written 



? (n) <k-D = 
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( A-16 ) 
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or 
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( A-17 ) 



where 



— (n) _ r -(n) -(n)-i 

w = L w . . .w J 

— xz xz n xz 

0 n 



( A-13 ) 



In the single channel case the backward weighting vector 
was equal to the forward weighting vector inverted in place. 
This arose because forward and backward reflection factors 
are equal in the single channel lattice. In the two channel 
lattice is not in general equal to , so there is no 

similar relationship between the forward and backward 
weighting vectors. 
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Using this notation 
substituting (A-12) and 



CA-^3) may be rewritten by 
(A-17), adjusted for proper order., 
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Comparing (A-19) with (A-12) and adjusting for order (n+1), 



; (n+ir 

-yy 



N <n+1) 

-yu 



w (n+1J l 

-uy 



„ <n+1) 

— UU 




-K 



( n+1 ) 




(A- 20) 



163 



This provides a recursive solution for the forward weight 
vectors . 

In illustration, e^ Q ^Ck) is defined by (A-l) as 
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(k) can be written, with a delay added, from (A-2) as 
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Substituting (A-21) and (A-22) into CA-20) yields 





r (i5i 

w 

-yy 










1 ' 
0 








1 

O i — 1 

.1 


n 


f 




(l) 

w 

[-yu J 










0 

_ 0_ 




T 

-K (1) 




r* 

0 

_ 0_ 








r w <^ 

-uy 










“ 0~ 








“ 0" 
















0 








0 








(1) 

w 

— uu 










1 

_ 0_ 








0 

_ 1_ 
















— 




- 















f 1 

i — 1 O 1 

1 I 








“0 “ 

k (1) 

*11 






1 


“0 - 
0 






L 


~i 

-k (1) 

*ii 






1 1 1 

1 1 

1 O O 1 

1 1 1 








1 ^ 

o o 

1 








0 

k (1) 

-21- 








0 

-k (1) 
- 21- 






0 

0 








"0 “ 

k (1) 

12 








1 

i 

0 O 1 

1 i 








“0 

-k (1) 

K 12 






1 1 

1 i — 1 O 

I 

1 






i 


1 ' 

o o 

1 








0 

k (1) 

L K 22_ 


1 


' 

1 


L 


1 

-k (1) 
- 2 2- 


1 



166 



This result can be substituted into (A~ll) 



e (1) (k) = 



e (1) (k) 

— 1 

e^ i; (k) 

u 



[y(k)y(k-l) i u(k)u(.k-l) ] 



y (k) -k^^y (k-1 ) -k^^uCk-l) 



_ 



-k^ 2 ^ y D +u ^ “ J<2 2 ^ U ^ 3<_ ^ 




-k 



( 1 ) 

22 



(A-24 ) 



Inspection of Figure 2.8 shows that (A-24) is indeed the 
equation of the first stage error signal. 

An analogous procedure can be allowed to develop the 
recursive solution for the backward weights. Substituting 
(A-12) and (A-17) into (A-4) yields 
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Comparing with (A-.17) and adjusting for order (n+1) reveals 
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With the weight vectors for all stages established, the 
reflection factors may be found in terms of these weights. 
Substituting (A-12) and (A-17) into (A-8) yields 
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The center matrix is dimensioned (2x2). Its four elements 
are each a submatrix of dimension C2 (n+2 )x2 (n+2 ) ) . Each 
submatrix is premultiplied by a 2 (n+2) element row vector 
and postmultiplied by a 2 (.n+2) element column vector. The 
result of these multiplications will be a scalar value for 
each of the four elements of A^ n ^* 

The order of these vector multiplications may be 
inverted yielding 
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Extracting A^-, and rewritting 
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is a (n+2)x(n+2) (cross) correlation array- 
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By defining a 2 (n+2 )x2 (n+2 ) matrix 
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Implementation is discussed in Chapter IV. B. 
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